#ifdef _WIN32
#  define NOMINMAX
#endif

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>

// includes, project
#include <helper_functions.h>
#include <helper_cuda.h>
#include "cuPrintf.cu"

// constants which are used in host and device code
#define             INV_SQRT_2      0.70710678118654752440f;
const unsigned int LOG_NUM_BANKS = 4;
const unsigned int NUM_BANKS = 16;

////////////////////////////////////////////////////////////////////////////////
// includes, kernels
#include "dwtHaar1D_kernel.cuh"

////////////////////////////////////////////////////////////////////////////////
// declaration, forward
void runTest(int argc, char **argv);
bool getLevels(unsigned int len, unsigned int *levels);

typedef long long int64;

struct PAIR {
	int a, b;
	PAIR(int a0, int b0) {
		a = min(a0, b0);
		b = max(a0, b0);
	}
};

bool operator<(const PAIR &x, const PAIR &y) {
	if (x.a == y.a)
		return x.b < y.b;
	else
		return x.a < y.a;
}
bool operator==(const PAIR &x, const PAIR &y) {
	return x.a == y.a && x.b == y.b;
}
struct hash_PAIR {
	size_t operator()(const PAIR &x) const {
		return (x.a << 8) ^ (x.b << 0);
	}
};
struct PII {
	int a, b;
	PII(int a0, int b0) {
		a = a0;
		b = b0;
	}
};
struct TRIPLE {
	int a, b, c;
	TRIPLE(int a0, int b0, int c0) {
		a = a0;
		b = b0;
		c = c0;
		if (a > b) {
			int temp = a;
			a = b;
			b = temp;
		}
		if (b > c) {
			int temp = c;
			c = b;
			b = temp;
		}
		if (a > b) {
			int temp = a;
			a = b;
			b = temp;
		}
	}
};
bool operator<(const TRIPLE &x, const TRIPLE &y) {
	if (x.a == y.a) {
		if (x.b == y.b)
			return x.c < y.c;
		else
			return x.b < y.b;
	} else
		return x.a < y.a;
}
bool operator==(const TRIPLE &x, const TRIPLE &y) {
	return x.a == y.a && x.b == y.b && x.c == y.c;
}
struct hash_TRIPLE {
	size_t operator()(const TRIPLE &x) const {
		return (x.a << 16) ^ (x.b << 8) ^ (x.c << 0);
	}
};

int n, m; // n = number of nodes, m = number of edges
int *n_device, *m_device;
int *deg; // degrees of individual nodes
int *deg_device;
PAIR *edges; // list of edges
PAIR *edges_device;

int **adj; // adj[x] - adjacency list of node x
int **adj_device;
PII **inc; // inc[x] - incidence list of node x: (y, edge id)
PII **inc_device;

int FLAG = 1;
int *flag_device;

int partition(int* a, int l, int r) {
	int pivot, i, j, t;
	pivot = a[l];
	i = l;
	j = r + 1;

	while (1) {
		do
			++i;
		while (a[i] <= pivot && i <= r);
		do
			--j;
		while (a[j] > pivot);
		if (i >= j)
			break;
		t = a[i];
		a[i] = a[j];
		a[j] = t;
	}
	t = a[l];
	a[l] = a[j];
	a[j] = t;
	return j;
}

void sort(int* x, int low, int high) {
	int j;

	if (low < high) {
		// divide and conquer
		j = partition(x, low, high);
		sort(x, low, j - 1);
		sort(x, j + 1, high);
	}
}

//TODO implement this partition kad je a == a, za b
int partition(PII* a, int l, int r) {
	PII pivot = a[l];
	int i = l;
	int j = r + 1;
	PII t = PII(0, 0);
	while (1) {
		do
			++i;
		while (a[i].a <= pivot.a && i <= r);
		do
			--j;
		while (a[j].a > pivot.a);
		if (i >= j)
			break;
		PII t = a[i];
		a[i] = a[j];
		a[j] = t;
	}
	t = a[l];
	a[l] = a[j];
	a[j] = t;
	return j;
}

void post_sort(PII* a, int n) {
	for (int i = 1; i < n; i++) {
		if (a[i - 1].a == a[i].a) {
			if (a[i - 1].b > a[i].b) {
				PII t = a[i];
				a[i] = a[i - 1];
				a[i - 1] = t;
			}
		}
	}
}

void sort(PII* x, int low, int high) {

	int j;

	if (low < high) {
		// divide and conquer
		j = partition(x, low, high);
		sort(x, low, j - 1);
		sort(x, j + 1, high);
	}
}

//TODO implement binary search
__device__ int binary_search(int* x, int n, int search) {
	int first = 0;
	int last = n;
	int middle = (first + last) / 2;

	while (first <= last) {
		if (x[middle] < search)
			first = middle + 1;
		else if (x[middle] == search) {
			return 1;
		} else
			last = middle - 1;

		middle = (first + last) / 2;
	}
	return 0;
}

int *adj_matrix; // compressed adjacency matrix
int *adj_matrix_device;

//returns boolean
__device__ int adjacent_matrix(int* adj_matrix, int x, int y, int* n) {
	int adj_chunk = 8 * sizeof(int);
	int k = adj_matrix[(x * (*n) + y) / adj_chunk]
			& (1 << ((x * (*n) + y) % adj_chunk));
	//cuPrintf("%d\n", k);
	return k;
}

int binary_search_host(int* x, int n, int search) {
	int first = 0;
	int last = n;
	int middle = (first + last) / 2;

	while (first <= last) {
		if (x[middle] < search)
			first = middle + 1;
		else if (x[middle] == search) {
			return 1;
		} else
			last = middle - 1;

		middle = (first + last) / 2;
	}
	return 0;
}

bool adjacent_list_host(int x, int y) {
	return 0; //binary_search_host(adj[x], adj[x] + deg[x], y);
}
const int adj_chunk = 8 * sizeof(int);

bool adjacent_matrix_host(int x, int y) {
	int k = adj_matrix[(x * n + y) / adj_chunk]
			& (1 << ((x * n + y) % adj_chunk));
	return k;
}

bool (*adjacent)(int, int);

int64 **orbit; // orbit[x][o] - how many times does node x participate in orbit o
int64 **orbit_device;
int64* orbit_array;

// triangles
int* tri_device;

//int *common;
//int *common_list;

//int64 *C4;
//int* neigh;

FILE *fin, *fout; // input and output files
int GS = 5;

int init(int argc, char *argv[]);

void count4_precomputing_host(int *tri, int64 *C4) {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac, frac_prev;

	// precompute triangles that span over edges
	printf("stage 1 - precomputing common nodes\n");
	//int *tri = (int*) calloc(m, sizeof(int));
	frac_prev = -1;
	for (int i = 0; i < m; i++) {
		frac = 100LL * i / m;
		/*if (frac != frac_prev) {
		 printf("%d%%\r", frac);
		 fflush(stdout);
		 frac_prev = frac;
		 }*/
		int x = edges[i].a, y = edges[i].b;
		for (int xi = 0, yi = 0; xi < deg[x] && yi < deg[y];) {
			if (adj[x][xi] == adj[y][yi]) {
				tri[i]++;
				xi++;
				yi++;
			} else if (adj[x][xi] < adj[y][yi]) {
				xi++;
			} else {
				yi++;
			}
		}
	}

	endTime = clock();
	printf("%.2f\n", (double) (endTime - startTime) / CLOCKS_PER_SEC);
	startTime = endTime;

	// count full graphlets
	printf("stage 2 - counting full graphlets\n");
//	int64 *C4 = (int64*) calloc(n, sizeof(int64));
	int *neigh = (int*) malloc(n * sizeof(int)), nn;
	frac_prev = -1;
	int k = 0;
	for (int x = 0; x < n; x++) {
		frac = 100LL * x / n;
		/*if (frac != frac_prev) {
		 printf("%d%%\r", frac);
		 fflush(stdout);
		 frac_prev = frac;
		 }*/
		for (int nx = 0; nx < deg[x]; nx++) {
			int y = adj[x][nx];
			if (y >= x)
				break;
			nn = 0;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = adj[y][ny];
				if (z >= y)
					break;
				k++;
				if (adjacent(x, z) == 0) {
					//fout << x << " " << z << "\n";
					continue;
				}
				neigh[nn++] = z;
			}
			for (int i = 0; i < nn; i++) {
				int z = neigh[i];
				for (int j = i + 1; j < nn; j++) {
					int zz = neigh[j];
					k++;
					if (adjacent(z, zz)) {
						C4[x]++;
						C4[y]++;
						C4[z]++;
						C4[zz]++;
					}
				}
			}
		}
	}

	endTime = clock();
	printf("%.2f\n", (double) (endTime - startTime) / CLOCKS_PER_SEC);
	startTime = endTime;
}

__global__ void count4_precomputing(int* m, int* n, int* deg, PAIR* edges,
		int** adj, int* tri_device) {

	cuPrintf("stage 1 - precomputing common nodes\n");

	int* tri = (int*) malloc((*m) * sizeof(int));
	int i = 0;
	//calloc
	for (; i < *m; i++) {
		tri_device[i] = 0;
	}

	for (i = 0; i < *m; i++) {
		int x = edges[i].a;
		int y = edges[i].b;

		int xi = 0, yi = 0;

		for (xi = 0, yi = 0; xi < deg[x] && yi < deg[y];) {
			if (adj[x][xi] == adj[y][yi]) {
				tri_device[i]++;
				xi++;
				yi++;
			} else if (adj[x][xi] < adj[y][yi]) {
				xi++;
			} else {
				yi++;
			}
		}

	}
}

//(m_device, n_device, deg_device, inc_device,
//		tri_device, flag_device, orbit_array, adj_matrix_device,
//		C4_device);
__global__ void count4_graphlets(int* m, int* n, int* deg, PII** inc,
		int* tri_device, int* flag, int64* orbit_array, int* adj_matrix,
		int64* C4) {
	int x = 2 * blockIdx.x;
	/*
	 // count full graphlets

	 //cuPrintf("stage 2 - counting full graphlets\n");

	 int64 *C4 = (int64*) malloc((*n) * sizeof(int64));
	 int *neigh = (int*) malloc((*n) * sizeof(int)), nn;

	 for (x = 0; x < *n; x++) {
	 C4[x] = 0;
	 neigh[x] = 0;
	 }

	 for (x = 0; x < *n; x++) {
	 for (int nx = 0; nx < deg[x]; nx++) {
	 int y = adj[x][nx];
	 if (y >= x)
	 break;
	 nn = 0;
	 for (int ny = 0; ny < deg[y]; ny++) {
	 int z = adj[y][ny];
	 if (z >= y)
	 break;
	 if (adjacent_matrix(adj_matrix, x, z, n) == 0) {
	 //cuPrintf("%d %d\n", x, z);
	 continue;
	 }
	 neigh[nn++] = z;
	 }
	 for (int i = 0; i < nn; i++) {
	 int z = neigh[i];
	 for (int j = i + 1; j < nn; j++) {
	 int zz = neigh[j];
	 if (adjacent_matrix(adj_matrix, z, zz, n)) {
	 C4[x]++;
	 C4[y]++;
	 C4[z]++;
	 C4[zz]++;
	 }
	 }
	 }
	 }
	 }
	 */
	//	cuPrintf("stage 3 - building systems of equations\n");
	int *common = (int*) malloc((*n) * sizeof(int));
	int *common_list = (int*) malloc((*n) * sizeof(int)), nc = 0;
	for (int x = 0; x < *n; x++) {
		common[x] = 0;
		common_list[x] = 0;
	}
	for (int x = 2 * blockIdx.x; x < 2 * blockIdx.x + 2; x++) {
		int64 f_12_14 = 0, f_10_13 = 0;
		int64 f_13_14 = 0, f_11_13 = 0;
		int64 f_7_11 = 0, f_5_8 = 0;
		int64 f_6_9 = 0, f_9_12 = 0, f_4_8 = 0, f_8_12 = 0;
		int64 f_14 = C4[x];
		orbit_array[x * 15 + 1] = 0;
		orbit_array[x * 15 + 2] = 0;
		orbit_array[x * 15 + 3] = 0;
		for (int i = 0; i < nc; i++)
			common[common_list[i]] = 0;
		nc = 0;

		orbit_array[x * 15] = deg[x];
		// x - middle node
		for (int nx1 = 0; nx1 < deg[x]; nx1++) {
			int y = inc[x][nx1].a, ey = inc[x][nx1].b;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = inc[y][ny].a, ez = inc[y][ny].b;
				if (adjacent_matrix(adj_matrix, x, z, n)) {
					if (z < y) {
						f_12_14 += tri_device[ez] - 1;
						f_10_13 += (deg[y] - 1 - tri_device[ez])
								+ (deg[z] - 1 - tri_device[ez]);
					}
				} else {
					if (common[z] == 0)
						common_list[nc++] = z;
					common[z]++;
				}
			}
			for (int nx2 = nx1 + 1; nx2 < deg[x]; nx2++) {
				int z = inc[x][nx2].a, ez = inc[x][nx2].b;
				if (adjacent_matrix(adj_matrix, y, z, n)) {
					orbit_array[x * 15 + 3]++;
					f_13_14 += (tri_device[ey] - 1) + (tri_device[ez] - 1);
					f_11_13 += (deg[x] - 1 - tri_device[ey])
							+ (deg[x] - 1 - tri_device[ez]);
				} else { // path
					orbit_array[x * 15 + 2]++;
					f_7_11 += (deg[x] - 1 - tri_device[ey] - 1)
							+ (deg[x] - 1 - tri_device[ez] - 1);
					f_5_8 += (deg[y] - 1 - tri_device[ey])
							+ (deg[z] - 1 - tri_device[ez]);
				}
			}
		}
		// x - side node
		for (int nx1 = 0; nx1 < deg[x]; nx1++) {
			int y = inc[x][nx1].a, ey = inc[x][nx1].b;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = inc[y][ny].a, ez = inc[y][ny].b;
				if (x == z)
					continue;
				if (adjacent_matrix(adj_matrix, x, z, n) == 0) {

					orbit_array[x * 15 + 1]++;
					f_6_9 += (deg[y] - 1 - tri_device[ey] - 1);
					f_9_12 += tri_device[ez];
					f_4_8 += (deg[z] - 1 - tri_device[ez]);
					f_8_12 += (common[z] - 1);
				}
			}
		}

		// solve system of equations
		orbit_array[x * 15 + 14] = (f_14);
		orbit_array[x * 15 + 13] = (f_13_14 - 6 * f_14) / 2;
		orbit_array[x * 15 + 12] = (f_12_14 - 3 * f_14);
		orbit_array[x * 15 + 11] = (f_11_13 - f_13_14 + 6 * f_14) / 2;
		orbit_array[x * 15 + 10] = (f_10_13 - f_13_14 + 6 * f_14);
		orbit_array[x * 15 + 9] = (f_9_12 - 2 * f_12_14 + 6 * f_14) / 2;
		orbit_array[x * 15 + 8] = (f_8_12 - 2 * f_12_14 + 6 * f_14) / 2;
		orbit_array[x * 15 + 7] = (f_13_14 + f_7_11 - f_11_13 - 6 * f_14) / 6;
		orbit_array[x * 15 + 6] = (2 * f_12_14 + f_6_9 - f_9_12 - 6 * f_14) / 2;
		orbit_array[x * 15 + 5] = (2 * f_12_14 + f_5_8 - f_8_12 - 6 * f_14);
		orbit_array[x * 15 + 4] = (2 * f_12_14 + f_4_8 - f_8_12 - 6 * f_14);
	}

}

__global__ void count5(int* m, int* n, int* deg, PAIR* edges, int** adj,
		PII** inc, int64 **orbit, int* flag, int64* orbit_array) {
	clock_t startTime, endTime;
	startTime = clock();
	clock_t startTime_all, endTime_all;
	startTime_all = startTime;
	int frac, frac_prev;

	int adj_chunk = 8 * sizeof(int);
	int size = ((*n) * (*n)) / adj_chunk + 1;
	int *adj_matrix = (int*) malloc(size * sizeof(int));
//calloc
	for (int i = 0; i < size; i++)
		adj_matrix[i] = 0;

	if (*flag == 1) {
		// compressed adjacency matrix
		for (int i = 0; i < *m; i++) {
			int a = edges[i].a, b = edges[i].b;
			adj_matrix[(a * (*n) + b) / adj_chunk] |= (1
					<< ((a * (*n) + b) % adj_chunk));
			adj_matrix[(b * (*n) + a) / adj_chunk] |= (1
					<< ((b * (*n) + a) % adj_chunk));
		}

	}

// precompute triangles that span over edges
	cuPrintf("stage 1 - precomputing common nodes\n");
	int *tri = (int*) malloc((*m) * sizeof(int));
	for (int i = 0; i < *m; i++) {
		tri[i] = 0;
	}
	frac_prev = -1;
	int k = 0;
	for (int i = 0; i < *m; i++) {
		frac = 100LL * i / (*m);
		if (frac != frac_prev) {
			//printf("%d%%\r",frac);
			//fflush(stdout);
			frac_prev = frac;
		}
		int x = edges[i].a, y = edges[i].b;
		for (int xi = 0, yi = 0; xi < deg[x] && yi < deg[y];) {
			if (adj[x][xi] == adj[y][yi]) {
				tri[i]++;
				xi++;
				yi++;
			} else if (adj[x][xi] < adj[y][yi]) {
				xi++;
			} else {
				yi++;
			}
		}
	}

	endTime = clock();
	cuPrintf("%.2f\n", (double) (endTime - startTime) / CLOCKS_PER_SEC);
	startTime = endTime;

// count full graphlets
	cuPrintf("stage 2 - counting full graphlets\n");
	int64 *C4 = (int64*) malloc((*n) * sizeof(int64));
	int *neigh = (int*) malloc((*n) * sizeof(int)), nn;
	for (int x = 0; x < *n; x++) {
		C4[x] = 0;
		neigh[x] = 0;
	}

	frac_prev = -1;
	for (int x = 0; x < *n; x++) {
		frac = 100LL * x / (*n);
		if (frac != frac_prev) {
			//printf("%d%%\r",frac);
			//fflush(stdout);
			frac_prev = frac;
		}
		for (int nx = 0; nx < deg[x]; nx++) {
			int y = adj[x][nx];
			if (y >= x)
				break;
			nn = 0;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = adj[y][ny];
				if (z >= y)
					break;
				k++;
				if (adjacent_matrix(adj_matrix, x, z, n) == 0) {
					//cuPrintf("%d %d\n", x, z);
					continue;
				}
				neigh[nn++] = z;
			}
			for (int i = 0; i < nn; i++) {
				int z = neigh[i];
				for (int j = i + 1; j < nn; j++) {
					int zz = neigh[j];
					k++;
					if (adjacent_matrix(adj_matrix, z, zz, n)) {
						C4[x]++;
						C4[y]++;
						C4[z]++;
						C4[zz]++;
					}
				}
			}
		}
	}

	endTime = clock();
	cuPrintf("%.2f\n", (double) (endTime - startTime) / CLOCKS_PER_SEC);
	startTime = endTime;

// set up a system of equations relating orbits for every node
	cuPrintf("stage 3 - building systems of equations\n");
	int *common = (int*) malloc((*n) * sizeof(int));
	int *common_list = (int*) malloc((*n) * sizeof(int)), nc = 0;
	for (int x = 0; x < *n; x++) {
		common[x] = 0;
		common_list[x] = 0;
	}

	frac_prev = -1;
	for (int x = 0; x < *n; x++) {
		frac = 100LL * x / (*n);
		if (frac != frac_prev) {
			frac_prev = frac;
		}

		int64 f_12_14 = 0, f_10_13 = 0;
		int64 f_13_14 = 0, f_11_13 = 0;
		int64 f_7_11 = 0, f_5_8 = 0;
		int64 f_6_9 = 0, f_9_12 = 0, f_4_8 = 0, f_8_12 = 0;
		int64 f_14 = C4[x];

		for (int i = 0; i < nc; i++)
			common[common_list[i]] = 0;
		nc = 0;

		orbit[x][0] = deg[x];
		// x - middle node
		for (int nx1 = 0; nx1 < deg[x]; nx1++) {
			int y = inc[x][nx1].a, ey = inc[x][nx1].b;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = inc[y][ny].a, ez = inc[y][ny].b;
				if (adjacent_matrix(adj_matrix, x, z, n)) {
					//if (adjacent(x, z)) { // triangle
					if (z < y) {

						f_12_14 += tri[ez] - 1;

						f_10_13 += (deg[y] - 1 - tri[ez])
								+ (deg[z] - 1 - tri[ez]);
					}
				} else {
					if (common[z] == 0)
						common_list[nc++] = z;
					common[z]++;
				}
			}
			for (int nx2 = nx1 + 1; nx2 < deg[x]; nx2++) {
				int z = inc[x][nx2].a, ez = inc[x][nx2].b;
				if (adjacent_matrix(adj_matrix, y, z, n)) {
					//if (adjacent(y, z)) { // triangle
					//		if (x == blockIdx.x)
					orbit[x][3]++;
					f_13_14 += (tri[ey] - 1) + (tri[ez] - 1);
					f_11_13 += (deg[x] - 1 - tri[ey]) + (deg[x] - 1 - tri[ez]);
				} else { // path
					orbit[x][2]++;
					f_7_11 += (deg[x] - 1 - tri[ey] - 1)
							+ (deg[x] - 1 - tri[ez] - 1);
					f_5_8 += (deg[y] - 1 - tri[ey]) + (deg[z] - 1 - tri[ez]);
				}
			}
		}
		// x - side node
		for (int nx1 = 0; nx1 < deg[x]; nx1++) {
			int y = inc[x][nx1].a, ey = inc[x][nx1].b;
			for (int ny = 0; ny < deg[y]; ny++) {
				int z = inc[y][ny].a, ez = inc[y][ny].b;
				if (x == z)
					continue;
				if (adjacent_matrix(adj_matrix, x, z, n) == 0) {

					orbit[x][1]++;
					f_6_9 += (deg[y] - 1 - tri[ey] - 1);
					f_9_12 += tri[ez];
					f_4_8 += (deg[z] - 1 - tri[ez]);
					f_8_12 += (common[z] - 1);
				}
			}
		}

		// solve system of equations
		orbit[x][14] = (f_14);
		orbit[x][13] = (f_13_14 - 6 * f_14) / 2;
		orbit[x][12] = (f_12_14 - 3 * f_14);
		orbit[x][11] = (f_11_13 - f_13_14 + 6 * f_14) / 2;
		orbit[x][10] = (f_10_13 - f_13_14 + 6 * f_14);
		orbit[x][9] = (f_9_12 - 2 * f_12_14 + 6 * f_14) / 2;
		orbit[x][8] = (f_8_12 - 2 * f_12_14 + 6 * f_14) / 2;
		orbit[x][7] = (f_13_14 + f_7_11 - f_11_13 - 6 * f_14) / 6;
		orbit[x][6] = (2 * f_12_14 + f_6_9 - f_9_12 - 6 * f_14) / 2;
		orbit[x][5] = (2 * f_12_14 + f_5_8 - f_8_12 - 6 * f_14);
		orbit[x][4] = (2 * f_12_14 + f_4_8 - f_8_12 - 6 * f_14);

		for (int j = 0; j < 15; j++) {
			//		cuPrintf("%llu ", orbit[x][j]);
			orbit_array[x * 15 + j] = orbit[x][j];
		}
		//	cuPrintf("\n");
	}

	endTime = clock();
	cuPrintf("%.2f\n", (double) (endTime - startTime) / CLOCKS_PER_SEC);

	endTime_all = endTime;
	cuPrintf("total = %.2f\n",
			(double) (endTime_all - startTime_all) / CLOCKS_PER_SEC);

}

void writeResults(int g = 5) {

}

void copyValuesToDevice() {
// initialize cuPrintf
	cudaPrintfInit();

	cudaMalloc((void**) &m_device, sizeof(int));
	cudaMemcpy(m_device, &m, sizeof(int), cudaMemcpyHostToDevice);

	cudaMalloc((void**) &n_device, sizeof(int));
	cudaMemcpy(n_device, &n, sizeof(int), cudaMemcpyHostToDevice);

	cudaMalloc((void**) &flag_device, sizeof(int));
	cudaMemcpy(flag_device, &FLAG, sizeof(int), cudaMemcpyHostToDevice);

	cudaMalloc((void**) &edges_device, m * sizeof(PAIR));
	cudaMemcpy(edges_device, edges, m * sizeof(PAIR), cudaMemcpyHostToDevice);

	cudaMalloc((void**) &deg_device, n * sizeof(int));
	cudaMemcpy(deg_device, deg, n * sizeof(int), cudaMemcpyHostToDevice);

	cudaMalloc((void**) &orbit_array, n * 15 * sizeof(int64));

	cudaMalloc((void***) &adj_device, n * sizeof(int*));
	int **adj_helper = (int **) malloc(n * sizeof(int *));
	cudaMalloc((void***) &adj_device, n * sizeof(int *));
	int i;
	for (i = 0; i < n; i++) {
		cudaMalloc((void**) &(adj_helper[i]), deg[i] * sizeof(int));
		cudaMemcpy(adj_helper[i], adj[i], deg[i] * sizeof(int),
				cudaMemcpyHostToDevice);
	}
	cudaMemcpy(adj_device, adj_helper, n * sizeof(int*),
			cudaMemcpyHostToDevice);

	free(adj_helper);

	PII **inc_helper = (PII **) malloc(n * sizeof(PII *));
	cudaMalloc((void***) &inc_device, n * sizeof(PII*));
	for (i = 0; i < n; i++) {
		cudaMalloc((void**) &(inc_helper[i]), deg[i] * sizeof(PII));
		cudaMemcpy(inc_helper[i], inc[i], deg[i] * sizeof(PII),
				cudaMemcpyHostToDevice);
	}
	cudaMemcpy(inc_device, inc_helper, n * sizeof(PII*),
			cudaMemcpyHostToDevice);
	free(inc_helper);

	int64** orbit_helper = (int64**) calloc(n, sizeof(int64*));
	cudaMalloc((void***) &orbit_device, n * sizeof(int64*));

	for (int i = 0; i < n; i++) {
		cudaMalloc((void**) &(orbit_helper[i]), 73 * sizeof(int64));
		cudaMemcpy(orbit_helper[i], orbit[i], 73 * sizeof(int64),
				cudaMemcpyHostToDevice);
	}
	cudaMemcpy(orbit_device, orbit_helper, n * sizeof(int64*),
			cudaMemcpyHostToDevice);
	free(orbit_helper);

	if (FLAG == 1) {
		int size = ((n) * (n)) / adj_chunk + 1;

		cudaMalloc((void**) &adj_matrix_device, size * sizeof(int));
		cudaMemcpy(adj_matrix_device, adj_matrix, size * sizeof(int),
				cudaMemcpyHostToDevice);
	}
}

void copyValuesFromDevice() {

	int mem_size = 15 * n * sizeof(int64);
	int64 orbit[15 * n];
	cudaMemcpy(orbit, orbit_array, mem_size, cudaMemcpyDeviceToHost);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < 15; j++) {
			if (j != 0)
				fprintf(fout, " ");
			fprintf(fout, "%llu", orbit[i * 15 + j]);
		}
		fprintf(fout, "\n");
	}

	fflush(fout);
	fclose(fout);

}

void freeValuesOnDevice() {

	cudaPrintfDisplay();
	cudaFree(m_device);
	cudaFree(n_device);
	cudaFree(edges_device);
	cudaFree(deg_device);
	cudaFree(adj_device);
	cudaFree(inc_device);
	cudaFree(orbit_device);

//todo free ostale promenljive

	cudaPrintfEnd();

}

void printOrbit() {
}

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
__global__
void saxpy(int n, float a, float *x, float *y) {
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if (i < n)
		y[i] = a * x[i] + y[i];
}

void effectiveBandwidth() {

	int N = 20 * (1 << 20);
	float *x, *y, *d_x, *d_y;
	x = (float*) malloc(N * sizeof(float));
	y = (float*) malloc(N * sizeof(float));

	cudaMalloc(&d_x, N * sizeof(float));
	cudaMalloc(&d_y, N * sizeof(float));

	for (int i = 0; i < N; i++) {
		x[i] = 1.0f;
		y[i] = 2.0f;
	}

	cudaEvent_t start, stop;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);

	cudaMemcpy(d_x, x, N * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_y, y, N * sizeof(float), cudaMemcpyHostToDevice);

	cudaEventRecord(start);

	// Perform SAXPY on 1M elements
	saxpy<<<(N + 511) / 512, 512>>>(N, 2.0f, d_x, d_y);

	cudaEventRecord(stop);

	cudaMemcpy(y, d_y, N * sizeof(float), cudaMemcpyDeviceToHost);

	cudaEventSynchronize(stop);
	float milliseconds = 0;
	cudaEventElapsedTime(&milliseconds, start, stop);

	float maxError = 0.0f;
	for (int i = 0; i < N; i++) {
		maxError = max((float) maxError, (float) abs(y[i] - 4.0f));
	}

	printf("Max error: %f\n", maxError);
	printf("Effective Bandwidth (GB/s): %f\n", N * 4 * 3 / milliseconds / 1e6);

}

int main(int argc, char **argv) {
	if (!init(argc, argv)) {
		printf("Stopping!");
		return 0;
	}
	copyValuesToDevice();
	cudaEvent_t start, stop;
	float milliseconds = 0;

	if (GS == 4) {
		int* tri = (int*) calloc(m, sizeof(int));
		int64* C4 = (int64*) calloc(n, sizeof(int64));
		count4_precomputing_host(tri, C4);
		cudaMalloc((void**) &tri_device, m * sizeof(int));
		cudaMemcpy(tri_device, tri, m * sizeof(int), cudaMemcpyHostToDevice);
		int64* C4_device;
		cudaMalloc((void**) &C4_device, n * sizeof(int64));
		cudaMemcpy(C4_device, C4, n * sizeof(int64), cudaMemcpyHostToDevice);

		cudaEventCreate(&start);
		cudaEventCreate(&stop);
		cudaEventRecord(start);
		count4_graphlets<<<n / 2, 1>>>(m_device, n_device, deg_device,
				inc_device, tri_device, flag_device, orbit_array,
				adj_matrix_device, C4_device);
		cudaEventRecord(stop);
		cudaEventSynchronize(stop);
		cudaEventElapsedTime(&milliseconds, start, stop);
	}
	if (GS == 5) {

		cudaEventCreate(&start);
		cudaEventCreate(&stop);
		cudaEventRecord(start);
		count5<<<1, 1>>>(m_device, n_device, deg_device, edges_device,
				adj_device, inc_device, orbit_device, flag_device, orbit_array);
		cudaEventRecord(stop);
		cudaEventSynchronize(stop);
		cudaEventElapsedTime(&milliseconds, start, stop);
	}
	copyValuesFromDevice();
	writeResults(GS);
	freeValuesOnDevice();
	printf("stage 3 - building systems of equations\n");
	printf("%f\n", milliseconds / 1000);
	printf("total: %f\n", milliseconds / 1000);
	return 0;
}

int init(int argc, char *argv[]) {
// open input, output files
	if (argc != 4) {
		fprintf(stderr, "Incorrect number of arguments.");
		fprintf(stderr,
				"Usage: orca [graphlet size: 4/5] [graph - input file] [graphlets - output file]");
		return 0;
	}
	sscanf(argv[1], "%d", &GS);
	if (GS != 4 && GS != 5) {
		fprintf(stderr, "Incorrect graphlet size %d. Should be 4 or 5.", GS);
		return 0;
	}
	char* mode = "r";
	fin = fopen(argv[2], mode);
	if (fin == NULL) {
		fprintf(stderr, "Can't open input file %s!\n", argv[2]);
		return 0;
	}

	char* writeMode = "w";
	fout = fopen(argv[3], writeMode);
	if (fout == NULL) {
		fprintf(stderr, "Failed to open file %s!\n", argv[3]);
		return 0;
	}

// read input graph
	fscanf(fin, "%d %d", &n, &m);
	int d_max = 0;
	edges = (PAIR*) malloc(m * sizeof(PAIR));
	deg = (int*) calloc(n, sizeof(int));

	for (int i = 0; i < m; i++) {
		int a, b;
		fscanf(fin, "%d %d", &a, &b);
		if (!(0 <= a && a < n) || !(0 <= b && b < n)) {
			fprintf(stderr, "Node ids should be between 0 and n-1.");
			return 0;
		}
		if (a == b) {
			fprintf(stderr, "Self loops (edge from x to x) are not allowed.");
			return 0;
		}
		deg[a]++;
		deg[b]++;
		edges[i] = PAIR(a, b);
	}
	for (int i = 0; i < n; i++)
		d_max = max(d_max, deg[i]);
	printf("nodes: %d\n", n);
	printf("edges: %d\n", m);
	printf("max degree: %d\n", d_max);
	fclose(fin);

//TODO FIX THIS
	/*if ((int)(set<PAIR>(edges,edges+m).size())!=m) {
	 cerr << "Input file contains duplicate undirected edges." << endl;
	 return 0;
	 }
	 */
// set up adjacency matrix if it's smaller than 100MB
	if ((int64) n * n < 100LL * 1024 * 1024 * 8) {
		//TODO FIX THIS
		FLAG = 1;
		adjacent = adjacent_matrix_host;
		adj_matrix = (int*) calloc((n * n) / adj_chunk + 1, sizeof(int));
		for (int i = 0; i < m; i++) {
			int a = edges[i].a, b = edges[i].b;
			adj_matrix[(a * n + b) / adj_chunk] |= (1
					<< ((a * n + b) % adj_chunk));
			adj_matrix[(b * n + a) / adj_chunk] |= (1
					<< ((b * n + a) % adj_chunk));
		}

	} else {
		FLAG = 2;
		//TODO FIX THIS
		adjacent = adjacent_list_host;
	}
// set up adjacency, incidence lists
	adj = (int**) malloc(n * sizeof(int*));
	for (int i = 0; i < n; i++)
		adj[i] = (int*) malloc(deg[i] * sizeof(int));
	inc = (PII**) malloc(n * sizeof(PII*));
	for (int i = 0; i < n; i++)
		inc[i] = (PII*) malloc(deg[i] * sizeof(PII));
	int *d = (int*) calloc(n, sizeof(int));
	for (int i = 0; i < m; i++) {
		int a = edges[i].a, b = edges[i].b;
		adj[a][d[a]] = b;
		adj[b][d[b]] = a;
		inc[a][d[a]] = PII(b, i);
		inc[b][d[b]] = PII(a, i);
		d[a]++;
		d[b]++;
	}
	for (int i = 0; i < n; i++) {
		sort(adj[i], 0, deg[i] - 1);
		sort(inc[i], 0, deg[i] - 1);
		post_sort(inc[i], deg[i]);
	}
// initialize orbit counts
	orbit = (int64**) malloc(n * sizeof(int64*));
	for (int i = 0; i < n; i++)
		orbit[i] = (int64*) calloc(73, sizeof(int64));
	printf("init done!\n");
	return 1;
}

////////////////////////////////////////////////////////////////////////////////
//! Perform the wavelet decomposition
////////////////////////////////////////////////////////////////////////////////
void runTest(int argc, char **argv) {
	bool bResult = false; // flag for final validation of the results

	char *s_fname = NULL, *r_gold_fname = NULL;
	char r_fname[256];
	const char usage[] =
			{
					"\nUsage:\n"
							"  dwtHaar1D --signal=<signal_file> --result=<result_file> --gold=<gold_file>\n\n"
							"  <signal_file> Input file containing the signal\n"
							"  <result_file> Output file storing the result of the wavelet decomposition\n"
							"  <gold_file>   Input file containing the reference result of the wavelet decomposition\n"
							"\nExample:\n"
							"  ./dwtHaar1D\n"
							"       --signal=signal.dat\n"
							"       --result=result.dat\n"
							"       --gold=regression.gold.dat\n" };

	printf("%s Starting...\n\n", argv[0]);

// use command-line specified CUDA device, otherwise use device with highest Gflops/s
	findCudaDevice(argc, (const char **) argv);

// file names, either specified as cmd line args or use default
	if (argc == 4) {
		char *tmp_sfname, *tmp_rfname, *tmp_goldfname;

		if ((getCmdLineArgumentString(argc, (const char **) argv, "signal",
				&tmp_sfname) != true)
				|| (getCmdLineArgumentString(argc, (const char **) argv,
						"result", &tmp_rfname) != true)
				|| (getCmdLineArgumentString(argc, (const char **) argv, "gold",
						&tmp_goldfname) != true)) {
			fprintf(stderr, "Invalid input syntax.\n%s", usage);
			exit(EXIT_FAILURE);
		}

		s_fname = sdkFindFilePath(tmp_sfname, argv[0]);
		r_gold_fname = sdkFindFilePath(tmp_goldfname, argv[0]);
		strcpy(r_fname, tmp_rfname);
	} else {
		s_fname = sdkFindFilePath("signal.dat", argv[0]);
		r_gold_fname = sdkFindFilePath("regression.gold.dat", argv[0]);
		strcpy(r_fname, "result.dat");
	}

	printf("source file    = \"%s\"\n", s_fname);
	printf("reference file = \"%s\"\n", r_fname);
	printf("gold file      = \"%s\"\n", r_gold_fname);

// read in signal
	unsigned int slength = 0;
	float *signal = NULL;

	if (s_fname == NULL) {
		fprintf(stderr, "Cannot find the file containing the signal.\n%s",
				usage);
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

	if (sdkReadFile(s_fname, &signal, &slength, false) == true) {
		printf("Reading signal from \"%s\"\n", s_fname);
	} else {
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

// get the number of decompositions necessary to perform a full decomposition
	unsigned int dlevels_complete = 0;

	if (true != getLevels(slength, &dlevels_complete)) {
		// error message
		fprintf(stderr, "Signal length not supported.\n");
		// cleanup and abort
		free(signal);
		exit(EXIT_FAILURE);
	}

// device in data
	float *d_idata = NULL;
// device out data
	float *d_odata = NULL;
// device approx_final data
	float *approx_final = NULL;
// The very final approximation coefficient has to be written to the output
// data, all others are reused as input data in the next global step and
// therefore have to be written to the input data again.
// The following flag indicates where to copy approx_final data
//   - 0 is input, 1 is output
	int approx_is_input;

// allocate device mem
	const unsigned int smem_size = sizeof(float) * slength;
	checkCudaErrors(cudaMalloc((void ** ) &d_idata, smem_size));
	checkCudaErrors(cudaMalloc((void ** ) &d_odata, smem_size));
	checkCudaErrors(cudaMalloc((void ** ) &approx_final, smem_size));
// copy input data to device
	checkCudaErrors(
			cudaMemcpy(d_idata, signal, smem_size, cudaMemcpyHostToDevice));

// total number of threads
// in the first decomposition step always one thread computes the average and
// detail signal for one pair of adjacent values
	unsigned int num_threads_total_left = slength / 2;
// decomposition levels performed in the current / next step
	unsigned int dlevels_step = dlevels_complete;

// 1D signal so the arrangement of elements is also 1D
	dim3 block_size;
	dim3 grid_size;

// number of decomposition levels left after one iteration on the device
	unsigned int dlevels_left = dlevels_complete;

// if less or equal 1k elements, then the data can be processed in one block,
// this avoids the Wait-For-Idle (WFI) on host side which is necessary if the
// computation is split accross multiple SM's if enough input data
	if (dlevels_complete <= 10) {
		// decomposition can be performed at once
		block_size.x = num_threads_total_left;
		approx_is_input = 0;
	} else {
		// 512 threads per block
		grid_size.x = (num_threads_total_left / 512);
		block_size.x = 512;

		// 512 threads corresponds to 10 decomposition steps
		dlevels_step = 10;
		dlevels_left -= 10;

		approx_is_input = 1;
	}

// Initialize d_odata to 0.0f
	initValue<<<grid_size, block_size>>>(d_odata, 0.0f);

// do until full decomposition is accomplished
	while (0 != num_threads_total_left) {
		// double the number of threads as bytes
		unsigned int mem_shared = (2 * block_size.x) * sizeof(float);
		// extra memory requirements to avoid bank conflicts
		mem_shared += ((2 * block_size.x) / NUM_BANKS) * sizeof(float);

		// run kernel
		dwtHaar1D<<<grid_size, block_size, mem_shared>>>(d_idata, d_odata,
				approx_final, dlevels_step, num_threads_total_left,
				block_size.x);

		// Copy approx_final to appropriate location
		if (approx_is_input) {
			checkCudaErrors(
					cudaMemcpy(d_idata, approx_final, grid_size.x * 4,
							cudaMemcpyDeviceToDevice));
		} else {
			checkCudaErrors(
					cudaMemcpy(d_odata, approx_final, grid_size.x * 4,
							cudaMemcpyDeviceToDevice));
		}

		// update level variables
		if (dlevels_left < 10) {
			// approx_final = d_odata;
			approx_is_input = 0;
		}

		// more global steps necessary
		dlevels_step = (dlevels_left > 10) ? dlevels_left - 10 : dlevels_left;
		dlevels_left -= 10;

		// after each step only half the threads are used any longer
		// therefore after 10 steps 2^10 less threads
		num_threads_total_left = num_threads_total_left >> 10;

		// update block and grid size
		grid_size.x =
				(num_threads_total_left / 512)
						+ (0 != (num_threads_total_left % 512)) ? 1 : 0;

		if (grid_size.x <= 1) {
			block_size.x = num_threads_total_left;
		}
	}

// get the result back from the server
// allocate mem for the result
	float *odata = (float *) malloc(smem_size);
	checkCudaErrors(
			cudaMemcpy(odata, d_odata, smem_size, cudaMemcpyDeviceToHost));

// post processing
// write file for regression test

	if (r_fname == NULL) {
		fprintf(stderr,
				"Cannot write the output file storing the result of the wavelet decomposition.\n%s",
				usage);
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

	if (sdkWriteFile(r_fname, odata, slength, 0.001f, false) == true) {
		printf("Writing result to \"%s\"\n", r_fname);
	} else {
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

// load the reference solution
	unsigned int len_reference = 0;
	float *reference = NULL;

	if (r_gold_fname == NULL) {
		fprintf(stderr,
				"Cannot read the file containing the reference result of the wavelet decomposition.\n%s",
				usage);
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

	if (sdkReadFile(r_gold_fname, &reference, &len_reference, false) == true) {
		printf("Reading reference result from \"%s\"\n", r_gold_fname);
	} else {
		cudaDeviceReset();
		exit(EXIT_FAILURE);
	}

	assert(slength == len_reference);

// compare the computed solution and the reference
	bResult = (bool) sdkCompareL2fe(reference, odata, slength, 0.001f);
	free(reference);

// free allocated host and device memory
	checkCudaErrors(cudaFree(d_odata));
	checkCudaErrors(cudaFree(d_idata));
	checkCudaErrors(cudaFree(approx_final));

	free(signal);
	free(odata);
	free(s_fname);
	free(r_gold_fname);

	cudaDeviceReset();
	printf(bResult ? "Test success!\n" : "Test failure!\n");
}

////////////////////////////////////////////////////////////////////////////////
//! Get number of decomposition levels to perform a full decomposition
//! Also check if the input signal size is suitable
//! @return  true if the number of decomposition levels could be determined
//!          and the signal length is supported by the implementation,
//!          otherwise false
//! @param   len  length of input signal
//! @param   levels  number of decomposition levels necessary to perform a full
//!           decomposition
////////////////////////////////////////////////////////////////////////////////
bool getLevels(unsigned int len, unsigned int *levels) {
	bool retval = false;

// currently signals up to a length of 2^20 supported
	for (unsigned int i = 0; i < 20; ++i) {
		if (len == (1 << i)) {
			*levels = i;
			retval = true;
			break;
		}
	}

	return retval;
}
